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I 

r 


i 


L_- 


n 


The  propagation  of  waves  in  both  fuUy  and  partially  ionized  gases,  both 
with  and  without  magnetic  fields,  has  been  treated  by  a  number  of  workers; 
e.  g.  Tanenbaum  and  Mintzer  in  1962  obtained  the  dispersion  relations  for  a 
linearized  and  spatially  uniform  gas  of  electrons,  positive  ions  and  neutrals. 
The  present  report  discusses  the  basic  formulation  and  mathematical  treatment 
of  wave  propagation  in  a  linearized  electron-ion  'neutral  gas,  with  static 
magnetic  field,  in  which  the  ambient  gas  parameters  have  an  arbitrary  varia¬ 
tion  with  the  vertical  coordinate  and  are  uniform  in  the  horizontal  direction. 

The  first  part  of  the  report  discusses  a  rather  standard  formulation  of  the 
general  problem,  via  the  Boltzmann  equation  and  the  Maxwell  equations.  By 
appropriate  momentum- space  averaging,  the  Boltzmann  equation  yields  motion, 
continuity  and  dynamic  adiabatic  state  equations.  These  are  then  combined  to 
yield  neutral  and  composite  plasma  equations  of  motion,  continuity  and  adia¬ 
batic  state  and  a  generalized  Ohm's  Law.  Steady  state  plane -wave  solutions 
are  suitable  in  the  horizontal  coordinates,  reducing  x,  y  and  t  dependence  to 
algebra  but  the  equations  must  remain  differential  in  the  vertical  coordinate  z. 
This  gives  rise  to  a  system  of  10  simultaneous,  highly  coupled  ordinary  first 
order  differential  equations  and  11  simultaneous  algebraic  equations. 

The  second  part  of  the  report  is  a  discussion  of  the  mathematical  solu¬ 
tions  of  this  coupled  algebraic  differential  equation  system,  which  is  equivalent 
to  the  system  of  equations  arising  in  analysis  of  coupled  linear  electrical  net¬ 
works.  Referring  to  both  the  mathematical  literature  of  differential  equations 
suid  the  modern  "state-space"  approach  to  automatic  control  systems,  various 


IV 


-2- 


purely  analytical  approaches  are  discussed  w>th  emphasis  on  their  deficiencies 
in  obtaining  practical  numerical  results  with  an  arbitrary  z-variation.  The 
Runge-Kutta  step-by-step  procedure  was  eventually  invoked  and  a  Fortran 
program  was  written  based  on  this  technique.  The  program  can  be  used  to 
obtain  accurate  numerical  solutions  to  many  problems  involving  waye  propaga¬ 
tion  in  a  linearized,  vertically  non-uniform  electron- ion-neutral  gas  without 
the  necessity  for  making  drastic  simplifying  assumptions  for  the  vertical  non¬ 
uniformity.  This  program  can  be  used  to  treat,  by  changing  input  parameter 
values,  such  diverse  problems  as  the  perturbing  effect  of  acoustic -gravity 
waves  on  ionospheric  electron  density,  electromagnetic  wave  propagation  in 
the  vertically  inhomogeneous  ionosphere,  MHD  waves  high  in  the  ionosphere, 
or  other  kinds  of  wave  propagation  in  plasra^a  media  with  a 
vertical  inhomogeneity.  Numerical  solutions  for  the  acoustic-gravity  wave 
plasma  interaction  problem  and  their  interpretation  will  be  presented  in  a  later 
report. 
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Dispersion  relations  for  a  homogeneous,  linearized  three  fluid  gas, 

consisting  of  electrons,  ions  and  neutrals,  with  a  static  magnetic  field,  have 

been  investigated  by  a  number  of  workers,  including  Tanenbaum  and  Mintzer^ 

2 

and  more  recently  by  Cronson  and  Clark  .  The  latter  focusses  attention  on 
the  specific  problem  of  the  ionosphere.  A  very  thorough  analytical  and 
numerical  studv  of  the  dispersion  relations  at  frequencies  below  a  few  cycles 
per  second,  using  the  latest  and  mbst  reliable  ionospheric  data  at  the  time  of 
writing  (June,  1964),  shows  that  of  the  four  permissible  modes  of  propagation, 
(neutral  acoustic,  plasma  acoustic  and  two  Alfven  modes),  the  only  ones  that 
are  not  prohibitive!}'  damped  out  by  collisions  between  charged  and  neutral 
particles  are  neutral  acoustic  waves.  Perturbation  of  the  plasma  gas  at 
appreciable  distances  from  the  exciting  source  at  low  frequencies  are  very 
likely  to  be  due  to  such  waves.  The  extreme  damping  of  other  plasma -induced 
propagation  modes  is  due  to  the  very  small  fractional  ionization  in  the  iono¬ 
sphere. 

Vertical  gradients  in  static  ionospheric  parameters,  e.  g.  neutral  gas 

density,  plasma  density  and  collision  frequency  are  in  some  cases  sufficiently 

steep  so  that  the  parameters  change  appreciably  within  a  wavelength  of  these 

low  frequency  acoustic  waves.  For  this  reason,  an  analysis  assuming  complete 

homogeneity  of  the  medium  might  exclude  some  important  effects.  In  particu- 
3  4  5  4 

lar,  Hines  ’  ’ '  and  Pitteway  have  analyzed  the  propagation  of  gravity  waves 
that  arise  in  a  neutral  gas  with  a  vertical  gradient  in  ambient  gas  pressure  and 
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density,  both  with  ’  and  without  consideration  of  viscosity  and  heat  conduc¬ 
tion.  An  analysis  not  accounting  explicitly  for  these  gradients  would  not  have 
predicted  the  ^-istence  of  such  waves. 

The  present  work  could  be  considered  as  a  generalization  of  Hines  1960 

3 

paper  to  include  the  interaction  of  the  gravity  waves  propagating  in  the  neutral 

gas  with  the  electron-ion  plasma  embedded  in  the  gas.  From  another  point  of 

view,  it  could  be  considered  as  a  generalizati  a  of  the  work  of  Tanenbaum  and 
1  2 

Mintzer  or  Cronson  and  Clark  to  include  the  effect  of  vertical  nonuniformity 
on  wave  propagation  in  a  composite  gas  consisting  of  electrons,  ions  and 
neutral  particles. 

The  present  paper  covers  the  formulation  of  the  general  problem  and  its 
mathematical  solution.  The  investigation  of  the  analytical  problem  that 
finally  led  to  the  use  of  a  numerical  technique  for  solution  is  described  in 
some  detail.  The  numerical  results  and  their  interpretation  will  be  treated 
in  a  later  report. 

2.  Formulation  of  the  General  Problem 

The  formulation  of  the  general  problem  (under  the  assumption  that  no 
ionization  processes  take  place  on  the  time  scale  of  the  phenomena  of  interest) 
begins  with  the  Boltzmann  equation  for  each  of  M  constituent  ideal  gases  and 
the  Maxwell  equations  for  electric  and  magnetic  fields.  These  are; 
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^  +  z  — 7^,  u'  ^  +  S  k  k 


(i) 


9t 


i=l  8x 


(i) 


1=1  ou  m 


'  V  9t  . 


;  k  =  1,  .  . .  M 


8h 

V  X  e  =  -  n 

9e 


(1) 

(2) 

(3) 


where 


til 

u;t)  =  number  density  of  particles  of  the  k —  species  at  a  point 
(r,  u)  in  phase -space  at  time  t, 

J  =  partial  time -rate  of- change  of  n^^  due  to  collision  effects^ 
c 

where  r  x^^^^  “epresents  position,  and 

u  =  u^^^^  the  velocity  vector,  where  "phase  space"  is  defined  as 

the  six-dimensional  space  of  position  and  velocity  (rather  than  position  and 
momentum). 

~  force  on  particles  of  k—  species, 
e  =  electric  field  vector  in  volts/meter 
h  =  magnetic  vector  in  amperes/ meter 
j  =  current  density  vector  in  amperes/meter^ 

-7 

fx  =  magnetic  permeability  of  free  space  =  47r  x  10  henry/meter 

1  -9 

€  =  diele^’tT’ic  coefficient  of  free  space  =  — —  x  10  farad/meter 

o  367r  ' 

Note  that  all  quantities  are  in  MK3  units. 

5 

in  the  standard  way,  following  Spitzer  ;  we  define  an  average  over  all 
velocity  space  of  a  general  function  G(u)  as; 
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G(u)  s  /  --  / G(u)n  (r,  uit)  d^u//  --  /  n  (r,u;t)d^u  (4) 

All  u- space  All  U“space 

and  note  that  the  number  density  of  k-type  par+icles  of  all  '/elocities  at  r  at 
time  t  is 


Nj^(r,t)  =  /  -  -  /  nj^(r,u;t)d\  (5) 

All  u- space 

By  setting  G(j)  =  1  in  (4),  we  are  led  to  an  equation  of  continuity  or  mass 

or  number  conservation  in  each  of  the  constituent  gases.  By  setting  G(u)  equal 

to  a  velocity  component,  we  are  led  to  equations  of  i.  otion  for  each  gas.  Two 

key  assumptions  in  these  developments  are: 

(a)  The  i—  component  of  the  average  force  F^^^is  independent  of  (i)  but 

may  depend  on  x^^^,  x^^^  and/or  u^^^  where  j  /  i.  Note  that  electric  and 

th 

gravitational  forces  are  independent  of  velocity  while  the  i—  component  of  the 
Lorentz  forces  ^  (u  x  h)  depend  only  on  velocity  components  other  than  u^^^ 


These  are  the  three  types  of  non-collision  forces  that  will  be  considered  in 

and  hence 


our  analysis.  Collision  effects  are  contained  m  tlie  term  ( 


are  not  included  in  the  force  F 

^  K 


(b)  The  number  density  N.  (r,t)  at  a  point  r  is  not  changed  by  collisions. 


1.  e. 


r  8Nj^(r;t) 

“at 


=  0 


c 


(«) 


where  the  subscript  c  denotes  "due  to  collisions". 

If  assumptions  (a)  and  (b)  are  invoked  in  (1)  and  the  averaging  process 
indicated  by  (4)  and  (5)  is  carried  out  with  G(u)  =  1,  the  result  is  a  set  of 
continuity  or  number-conservation  equations. 
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V  •  (N,  u  )  +  =  0;  k  =  1,  .  .  .  ,  M  17) 

k*k  9t 

tilt 

where  is  the  average  velocity  of  particles  of  the  k —  species  or  equivalently 

®‘’k 

V  •  (P^u^>  +  -g^  =  0;  k  =  1 . M  (V) 

where  p  h  N  ,  the  mass-density  of  k-type  particles. 

iC  K  iC 

Application  of  the  same  averaging  process  with  the  assumptions  (a)  and 
(b)  to  the  case  where 

G(u)  =  u^^^  i  =  1,  2,  3  (8) 

th 

results  in  the  vector  equation  of  motion  for  the  k —  particle  species,  i.  e. 


^k  Dt  ^k  (  8t  -^k  — k 


r^  +  -^(e  +  p  u  xh) 
-  mk  -*  o-k 


+  p  V0  +  Z'  p  V  .(u  -  u.) 

k  ^g  ^k  kj-k  -3 


where  the  forces  accounted  for  are  electric,  magnetic,  and  gravitational  field 

forces,  the  latter  being  proportional  to  the  gradient  of  a  gravitational  potential 

th 

function  <f)^.  q^^  's  the  charge  of  a  particle  of  the  k —  species;  the  collision 
forces  are  assumed  to  be  of  the  smiple  '  irictional”  type,  as  obtained  from 
simple  kinetic  theory  arguments,  i.e.  proportional  to  the  velocity  difference 


between  colliding  particles.  The  collision  frequencies  are  given  by  v  and 

k] 

(k) 

the  stress  tensor  ijf  is  defined  as  a  matrix  whose  elements  are: 

Assumption  of  a  Maxweil-Boltzmann  distribution  of  velocities  diagonalizes  the 


stress-tensor  and  by  very  elementary  kinetic  the  ory  arguments,  we  obtain 
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—  ,.,,2 


where  v,  =  L 
speed 


r  U)  U)T 

L  u  -  u  =  mean  square  deviation  of  k-type  particle 

_  2  _  K  R 


p  =  scalar  pressure  of  k —  gas  species. 

JV 

Auxiliary  to  eq.  (9)  is  the  statement  that  tlie  net  momentum  exchange  in 
collisions  between  any  two  types  of  particles  must  be  zero,  required  by 
Newton's  third  law.  Thus; 

Vkj<“k-“j>'- . 

as  a  consequence  of  which 
V,  .  p.  m.N. 

.  (13 

'’jk  •’k  ”k”k 

We  now  perform  the  standard  linearization  procedure  on  che  equations 
(2),  (3)  (7)  or  (7)'  and  (9);  i.r.  we  assume  that  the  variables  N.  or  p  ,  p 

R  R  R 

and  the  components  of  h  each  consist  of  a  large  zero-order  part,  denoted  with 

subscript  0,  plus  a  small  first  order  perturbation  part,  denoted  with  subscript 

1,  while  the  components  of  u^  and  e  are  assumed  to  have  vanishing  first  order 

parts,  i.  e.  to  be  of  "small  perturbation"  magnitude.  We  then  neglect  all  terms 

of  second  order,  e.  g.  the  term  u  •  Vu  in  (9)  or  the  term  V  •  (p,  ,u  )  in  (7)'. 

^  R  ^  R  R  X^R 

We  then  assume  that  the  pressure  and  density  of  each  gas  are  related  by  an 
adiabatic  state  equation  whose  general  and  linearized  forms  are: 
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2  ^k1  - 

~  ■  ‘=sko  -at-  =  ■  “k  '  I'^Pko  ■  ^sko'^Pko!  (l^^^rlzed) 
where  y  is  the  ratio  of  specific  heat  at  constant  pressure  to  that  at  constant 


volume,  assumed  to  be  the  same  for  all  gases,  and  c 


2  YPko  . 


is  the  zero 


order  sound  velocity  in  the  k —  gas  species.  Eqn.  (14)  is  based  on  the 
assumption  that  no  heat  flows  in  or  out  during  the  short  periods  associated 
with  the  perturbation  of  the  gas,  although  the  ambient  or  "zero  order”  state 
of  the  gas  is  by  no  means  necessarily  adiabatic. 

Finally,  we  write  the  set  of  linearized  equations  required  to  describe  the 


system,  dropping  the  lines  over  the  u^'s  which  indicate  averaging,  as  follows; 


®Sk  Pko^ 

1-1^:  =  -  ^P,.,  +  (e  +  «  u,  xh  )  +  p,  ,  VA 

'^ko  at  *^kl  m^^  -  ^o  ^kl  ^g 


•".='Pko’'kA'^j*'‘‘'^ . 


Eq.(9) 


‘"kolik*  ^  ° 


m-k; 


2  ^Pkl 
sko  at 


Eq.(7)' 


“k'  t'^>^.o-‘=sko  VpJ:k  =  l,...,M  Eq.(14)' 


IV;  V  X  e  =  -  M 


o  at 


Eq.(2) 


V;  Vxhj=j  + 

where  the  current  density  j  is  defined  by 


Eq.(3) 


j  =  2  j 
-  k=l-f' 


where  1  ^  N^^q^a^  = 


"ko'lk'ik 


Note  that  the  condition  expressed  by  eq.  (13)  must  hold. 


-8- 


) 

I)' 

)' 

) 


I) 


Each  of  the  vector  equations  I-k,  IV  and  V  has  3  component  equations 
(which  we  will  denote  by  eq.  I-k-1,  I-k-2^  I-k-3j  IV-1,  etc,).  Thus  the 
system  I-k  through  V  contains  a  total  of  (3M  +  M  +  M  +  3  +  3)  =  (5M  +  6) 
equations.  The  zero  order  parameters  are  assumed  to  be  known.  The  un¬ 
knowns  are  the  perturbed  or  first  order  parameters,  i.  e.  all  components  of 
^'s,  constitucing  a  total  of  3M,  all  of  che  which 

constitutes  a  total  of  M  parameters,  and  ihe  components  of  e  and  h,  which 
adds  another  six.  The  total  number  of  unknowns,  then,  is  (5M+6),  which  is 
the  same  as  the  number  of  parameters.  The  components  of  j  are  not  addi- 

tional  unknowns,  because  i  is  a  linear  combination  of  the  u  's. 

-k 

We  now  note  that  each  gas  was  assumed  to  obey  an  ideal  gas  law 


Is 

lo 


tf.o 


JLO 


=  kT 


lo 


(16) 


where  k  is  Boltzmann's  constant  and  T^^  the  zero-order  absolute  temperature 
til. 

of  the  I —  constituent  gas. 


3.  The  Neutral-Electron-Ion  Gas 

We  will  now  specialize  to  the  three-fluid  (neutral,  positive  ion,  electron) 

gas,  the  homogeneous  case  of  which  was  treated  by  Tanenbaum  and  Mintzer^ 

2 

and  by  Cronson  and  Clark  .  In  this  case: 

k  =  1;  neutral  gas;  subscript  n  replaces  subscript  1  on  all  quantities; 


k  =  2;  ion  gas,  subscript  i  replaces  subscript  2  on  all  quantities;  q,,  =  +  eZ, 

where  e  is  the  electron  charge’!'  and  Z  the  degree  of  ionization. 

Not  to  be  confused  with  e,  the  electric  field,  or  its  components,  although 
the  same  letter  is  used  to  represent  both  quantities. 
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k  =  3;  electron  gas;  subscript  e  replaces  subscript  3;  q  =  -  e. 

O 

We  now  invoke  the  following  definitions  and  assumptions,  most  of  which 

2 

were  used  by  Cronson  and  Clark  , 

Definitions: 

p  s  partial  pressure  of  plasma  (electron  plus  ion  gas)  =  p.  +  p 
p  1  e 

Pp  =  mass  (density  of  plasma  =  ^e^e 

p.  u.  +  p  u 
1  1  io-»i  eo-^e 

u  =  average  plasma  velocity  =  - - - 


lO 


j  =  eZN.  u.  —  eN  u 
d.  10*1  eo^e 


<m 


m. 

1 


Assumptions: 


P.  P 

- Z  =  N.  Z  =  N  = -  ;  electrical  neutrality  (in  zero-order) 

m.  10  eo  m  '' 

i  e 


(17) 


from  which  it  follows,  through  the  definition  of  j  above,  that 

ep. 


eo 


(u,  -  u  ) 


m  ^1  j.e 

e 


3  = 

Zm 

0 

- - «  m.,  or  ^  «  1;  (Electron  mass  is  extremely  small  compared  to  ion 

i 

mass  and  Z  is  a  number  of  order  unity) 


(18) 


(19) 


from  which  it  follows,  through  the  above  definition  of  u  ,  that 

-►P 

N  m 
,  eo  e 
u  =  u.  +  - —  - u 


.p  ^i  N.  m.  *e  u.  +  ^  u 
lo  1  ^1 


1  + 


N  r 
eo 

N.  m. 

lO  1 


^  j.  t  ~  u.  +  I  u 
1+1  *1  ^  *e 


(20) 


From  the  ideal  gas  law. 


T  p  N. 
eo  _  eo  lo 

T.  ""  p.  N 
lo  lo  eo 


=  1 


(21) 
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i.  e.  and  the  zero-order  electron  and  ion  temperature  are  assumed  equal, 
or  equivalently,  thermal  equilibrium  between  electrons  and  ions  is  assumed. 

We  now  note  that  (see  eqn.  11) 


N  m 
eo  e 


Ip 


1  i  N.  m 

^  lO 


l“e  [klisl 


■ 

2 


(22) 


and  assume  that  the  electrons,  being  many  times  lighter  than  the  ions,  attain 
much  higher  rms  fluctuation  velocities.  This  assumption  results  in  the  state¬ 
ment  that 


p  »  I  p. 
e  1 

/Iso,  from  (17)  and  (21) 


(23) 


eo 

Pio 


N 

eo 

N. 

10 


=  Z 


Note  that,  from  (24) 
2 


seo 


sio 


P  P 

eo  ^io  _  _Z_  j 

P.  P  ~  I 

10  eo  ^ 


We  now  define  a  "plasma  sound  velocity"  by 


2  2 
YP  c  .  p .  -P  c  p 
po  sio  10  seo  *^00 


spo 


po 


p.  +  p 
lo  eo 


^’PjQ  ^  Peo> 

<Plo  Peo> 


(24) 


(25) 


^sio  «  c  ^(1  +  Z) 

sio 


(l+l) 


and  note  that 


(26) 


spo 


sno 


^  . 


P 


no 


a 


(27) 


where 
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*^no 


=  fractional  ionization. 


By  solving  (18)  and  (20)  simultaneously  for  u.  and  u  ,  with  the  aid  of  (19). 

X  6 


we  obtain 


u.  a 
.^1 


u  ~ 


_*e 


u  + 

-P 

u  - 


3 


3 


(28) 

(29) 


where  ri  =  -  . 

m 

e 

We  now  write  the  equations  I-k,  Il-k,  Ill-k  for  the  case  under  considera¬ 


tion 
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e 
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TTI  1  ^  ‘  U  )  +  1  ~  ® 

II’ -1  no  n  \9t /  nl 
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1  1 
e  e 


1 
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OjJl 

i 
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m'-i 

-2 

-3 
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^Inl-  sno 

'  i  i 


(IfV. 


u  [Vp  -  c  Vp  ] 
-.^n  •  no  sno  '^no 
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1 
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1 

e 


1 

e 


1 
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The  following  operations  are  now  performed; 

(A)  Add  I' -2  to  I' -3 

(B)  Subtract  (I' -3  multiplied  by— ^)  from  (I' -2  multiplied  by^^^) 

m  m. 

e  1 

(C)  Add  II' -2  to  II' -3 

(D)  Subtract  (II' -3  multiplied  by  — ^)  from  (II' -2  multiplied  by  — ) 

e  "^i 


(E)  Add  III' -2  to  nr -3 
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(F)  Subtract  (III' 3  multiplied  by - )  from  IIP  2  multiplied  by  - ) 

m  m. 

e  1 

Performance  of  these  operations  and  use  of  the  definitions  of  p  ,  p  ,  u 

P  P  -P 

and  j  given  above,  eq.  (17),  together  with  eqs.  (13)  and  eqs.  (17).  (19).  (21), 

■* 

(23),  (25),  (26)  (28)  and  (29).  results  in  the  following  set  of  equations  describ¬ 


ing  the  sy&cem. 


T"  V, 

I  -n,  p 


9u 

.-n 

no  9t 


^P  H  P  1  +  ctp  (v .  +  £  V  )  (u  -  u  )  + 

nl  nl  ^g  '^no  in  ^  en  Un  ^.p 

(Motion  equation  for  neutral  gas) 


V  TV.  V 

en  .inN 

y 


9u 

'"■P-  Ppo  if  ° 


Ip  (v.  +£v  )(u  -u) 

^pl  ^pl  i  ,aO  '^po  m  ^  en  ^p  >n 


-V.  -  V  V 


(where  B  has  replaced  p.  h  ) 
^  o  o,*o 


‘"-'’if  = 


2,^, 

e  N 


^Vp  +  - 

1  m 


eo 


+  T1P  £(v  -v.)(u-u)4-(v  -P|v.  +v.)i-i  (gp  , ) 

•'^po^  en  m  ^n  ^p  en  ^  in  ei  1  z  ®*^cl 


e^N 


.e^N 


+  - -  u  XB  -  - - ^  (j  xB  ) 


where  p^  Pil"  Pgi^ 

(generalized  Ohm's  Law  for  plasma) 

9]P  . 

II''-n  V  •  (p  u  )  +  .P  =  0  (continuity  for  neutral  gas) 
noun  ot 


9p  1 

II" -p:  V  •  (p  u  )  +  — r“  =  0  (continuity  for  plasma) 
po-*p  9t 

Op 


II" -0:  V  •  j  + 


cl 


9t 


=  0  (electrical  equation  of  continuity  or  charge  conservation) 
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where  =  perturbed  or  first  order  charge  density, 

not  necessarily  zero 

9p„i  9  o 

IIl"-n:  — TT —  -  c  -r- —  =  -  u  •  [Vp  -  c  Vp  ] 

ot  sno  8t  ^n  no  sno  no 

(Adiabatic  state  for  neutral  gas) 

in"-p:  -  c  ^  =  -  u  •  [  Vp  -  c  ^Vp  ] 

8t  spo  at  -^P  po  spo  po^ 

_  r  _  2  .  Z  *^spc  ^^cl 

^  “nppod  +  Z)  ^po  ^spo  Ppo^  (1+Z)  ^r\  at 

(Adiabatic  state  for  plasma) 


TTTii  n  1  .  ^^spo  ^^pl  Z  r  „  1 

III  -0:  -3^  +  Up  •  [  Vpp_^  -  Vppj 

2  . 

_  i _  •  r  V7  -  ^\7  14.  '^i^po  *^cl 

(]+Z)t)|p  ^no  ^spo  ^po  |ti  at 


jv".  V  X  e  =  -  M 


o  at 


V":  Vxh^=j  +  €^  -5= 


Note  that  there  are  other  equally  satisfactory  ways  to  describe  the  three - 

fluid  (electron-ion-neutral)  gas,  such  as  the  retention  of  the  original  equations 

r-1,2,3,  II' -1,2,3,  III' -1,2,3,  wherein  the  unknowns  are  u  .  p  .  - 

n,  1,  e,  n,  1,  e. 

etc.  ,  instead  of  u  ,  p  ,  j,  as  in  the  present  formulation.  Tanenbaum  and 

n, p  n, p 

1  2 
Mintzer  use  the  (n,  i,  e. )  form  of  the  equations,  while  Cronson  and  Clark  , 

g 

following  Spitzer  ,  use  the  ^n,  p)  form.  The  latter  approach  conveniently  lumps 
the  electron  and  ion  gas  equations  into  a  single  plasma  equation,  and  the  current 
density  j  appears  in  the  Maxwell  equations  as  a  single  unknown  vector  instead  of 
a  linear  combination  of  the  vectors  u.  and  u  .  However,  there  is  no  real 
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difference  in  principle  or  in  mathemat.ic?\l  difficulty  between  these  two  methods 
of  approach.  Both  leave  the  same  number  of  equations  and  unknowns,  but  the 
equations  appear  different  in  form  and  the  unknowns  are  changed. 

Still  another  method,  that  used  by  Cowling^,  lumps  the  three  gases 
together  into  a  composite  neutral-plus -plasma  gas,  defining  variables  for  that 
composite  gas.  The  approach  illuminates  different  physical  parameters  than 
the  present  approach,  but  does  not  change  the  degree  of  difficulty  in  solving  the 
equations.  Moreover,  solutions  obtau^ed  with  one  of  these  approaches  can  al¬ 
ways  be  used  to  obtain  the  solutions  that  would  have  been  obtained  by  another, 
since  the  former  sets  of  solutions  are  linear  combinations  of  the  latter. 

4.  The  Vertically  Non-uniform  Neutral-Electron-lon  Gas 

We  now  further  specialize  to  the  case  where  the  ambient  three -fluid  gas 
has  parameters  that  are  constant  in  time  and  horizontally  uniform,  but  vary 
with  vertical  position.  To  treat  this  case,  we  use  a  right-handed  cartesian 
coordinate  system  (x,  y,  z)  and  look  for  plane  wave  solutions  in  the  time  t  and 
the  horizontal  coordinates  x  and  y,  reducing  the  problem  to  a  system  of  ordinai'y 
differential  equations  in  the  vertical  coordinate  z.  Thus  each  first  order  (un¬ 
known)  parameter,  designated  generically  by  u,  has  a  complex  Fourier  solution 
of  the  form 

-i(wt  -  k  X  -  k  y) 

u(x,  y,  z,t)  =  U(k  ,k  z,a))e  ^  ^  ,  (30) 

X  y 

(i,  e.  the  wave  solution  corresponding  to  a  given  lower  case  letter  is  denoted 

by  its  corresponding  capital  letter,  except  for  the  solution  corresponding  to  p, 

which  is  denoted  by  01)  where  k  and  k  are  in  general  complex  while  the  angular 

X  y 

frequency  u>  is  real. 


15- 


Wo  define  a  horizontal  wave  vector 


kj^by 


k,  =  i  k  +  i  k 
X  y 


(31) 


where  are  unit  vectors  in  the  x,  y,  z  directions  resulting  in  the 

following  forms  for  eqs.  I"  through  V". 
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111“' -n]  -  i«P  +  c  i()(ft  =  -  U 

nl  sno  nl  nz 
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IX 
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It  is  convenient  to  normalize  the  quantities  in  the  equations  I“'-n  through 
V'  such  that  each  quantity  in  the  equations  is  dimensionless.  To  this  end  we 
denote  each  normalized  parameter  by  its  symbol  with  a  karet  above  it  and  in¬ 
voke  the  following  aefinitions: 


en,  in,  ei 


en,  in,  ei 

CJ 


sko 


sko 


;  k  =  n,  p,  i,  e 


c  =  light  velocity  in  vacuo  = 


^  ti  e 
o  o 


u  =  'HP  H  =  electron  cvclotron  frequency 
ce  o  o 
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K  .  = 


sko 


ck  c  ,  3z 

sko 
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k  =  free  space  electromagnetic  wave  number  =  —  =  co*/S  e 
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X  k 


CO  =  electron  plasma  resonant  frequenc/  = 
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/  we  ju  V 

B,  =„  H 
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V,  =  (  -r  -  V, 


i-j  r  =  unit  vector  in  direction  of  H 
-B  H  ^o 

'  j.o' 
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cl  T)^  p 

DO 


We  also  assume  that; 

V<|)  =  g  =  -gi  =  -  gwci  (32 

g  ^  «•  z  ^  z 

m  s 

where  g  =  acceleration  to  gravity  =  9.  8 - ;r  and  g  =  — ^  . 

2  wc 

sec 

and  conclude  from  (23)  that 

P.  =  I  P-  -  P  ~  -  P  (33 

i  1  e  e 

i.  e.  the  parameter  p^,  which  enters  into  later  discussions,  is  approximately 
equated  to  the  negative  of  electron  pressure. 

In  terms  of  the  normalized  parameters  and  with  the  aid  of  the  assumption 
(32),  we  have 
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4.  1.  The  Two-Dimensional  Caue 

Consider  the  case  where  the  source  of  the  wave  disturbance  under  study 
is  a  uniform  infinitely  long  line-source  along  the  y-axis  and  the  ambient  gas 
parameters  are  independent  of  x  and  y.  With  this  kind  of  symmetry,  the  un¬ 
knowns  cannot  be  functions  of  y.  We  can  therefore  set  k  =  0  and  the  horizontal 
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propagation  vector  becomes  We  now  resolve  the  vector  equations 

A  A  A 

I'".  IV'"  and  V"  into  component  equations,  assigning  the  aopropriate 
symbols  to  denote  component  equations  x,  y  and  z. 
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4.  1.  1.  Matrix  Formulation 


Eqns.  I-n-x  through  V-z  constitute  a  set  of  21  coupled  equations  in 
21  unknowns,  of  which  11  are  purely  algebraic  and  10  are  first  order 
differential.  To  solve  them  by  matrix  methods,  we  first  express  them  in 


the  form 

,  =  A'(z)  u(z)  +  B'(z)  v(z) 

dz 

(34) 

C  =  C'(z)  u(z)  +  D'(z)  v(z) 

(35) 

where  A'(z).  B'(z),  C'(z)  and  D'(z)  are  matrices  consisting  of  known  elements 
composed  of  certain  coefficients  of  our  equation  system,  while  u(z)  and  v(z) 
are  vectors  whose  elements  are  the  unknown  parameters.  The  matrices  A’, 
B’,  C  and  D'  and  the  vectors  of  u  and  v  are  shown  in  Tables  1  and  2. 

To  cast  Eqns.  (34)  and  (35)  into  the  desired  form,  we  perform  the 
following  matrix  operations: 

(D')’^C'u  +  D'v)  =  0  =  (D')'^  C'u  +  V  ; 


or 


v  =  -  (D')’^  C'u 


du 

dz 


=  A'u  +  B'  [-  {(D')’^  C'}u] 


[A'  -  B'(D')’^C']u 


or  more  concisely 


du 

dz 


=  Au 


(36) 

(37) 

(37)' 


where 
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Note  that  th“  operation  indicated  in  matrix  notation  by  Eq.  (36)  is  that  of 
solving  the  set  of  11  simultaneous  eq;iai;ions  for  the  elements  of  v,  considered 
unknown,  in  tetii^a  of  those  of  u,  considered  known.  The  next  step  represents 
the  subsequent  operation  of  substituting  these  solutions  of  (34)  iixto  (35)  in 
order  to  obtain  a  system  of  coupled  first  order  d  rferential  equations  of  the 
form  (37). 


4. 1..  2.  Analytical  Methods  of  Solution 


The  theory  of  systems  of  the  class  (37)  is  discussed  in  both  the  purely 

(8  9  10  11  12) 

mathematical  literature  of  differential  equations'  •  •  ‘  ’  and  in  the 

literature  of  mod'^-'n  generalized  systems  and  automatic  control  theory 
(13,  14,  15,  16)  particular,  the  system  (37)  with  constant  A  elements  is 
the  foundation  of  the  "state-space*'  approo.ch  to  modern  automatic  control 
theory^ The  literature  was  reviewed  in  an  effort  to  find  techniques 
for  solution  of  the  syst'^m  (37)  that  apply  to  matrices  A(z)  with  z-dependent 
elements.  If  these  elements  are  independent  of  z,  or,  even  if  z-depsndent, 
if  they  have  the  property  of  seK-comrx.utativity,  i.e. 

A(z  )  A(2.,)  =  A(z  )  A(z  )  (38 

1  z  1 

then  solutions  can  be  obtained  in  closed  form^^^'  the  problem  degenerates 
into  that  of  finding  eigenvalues  of  the  matrix  /  A(z)dz.  If  A(z)  does  not  have 
the  property  (38),  then  each  of  the  elements  r.- -y  be  separable  into  a  major 
part  that  fulfills  (38)  and  a  "small  perturbation"  part  that  does  not.  The 
system  is  then  solved  with  the  assumption  that  the  entire  matrix  is  the  self- 


commutative  part,  and  this  solution  is  used  as  the  zero-order  tferm  of  a 


perturbation  solution  for  the  more  realistic  problem  involving  the  non  self- 
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commutative  matrix.  The  difficulty  with  this  approach  is  that  it  is  very 
difficult  to  find  a  meaningful  separation  of  A(z)  into  self-commutative  and  xion 
self-commutative  parts. 


Another  possible  analytical  technique  applies  to  periodic  elements  of 
(8  9  10  11) 

A(z)'  "  '  ■  .If  the  z -range  is  divided  into  strata,  then  in  each  stratum 

the  elements  A(z)  could  be  considered  periodic  with  a  period  equal  to  the 

(9) 

size  of  the  stratum.  General  theory  exists  on  the  periodic  coefficient  case  , 
but  its  implementation  requires  knowledge  of  the  "fundamental  solutions"  of 
the  system,  which  is  not  easily  attainable,  or  what  amounts  to  a  perturbation 
around  the  constant  coefficient  case,  which  may  be  very  maccurate  unless 
the  coefficients  are  nearly  constant  or  unless  many  terms  of  a  Fourier  series 


solution  are  invoked. 


Another  possibility  that  was  studied  is  the  straightforward  diagonalization 
of  A(z)  which  effectively  uncouples  the  differential  equations  (37)  and  trans¬ 
forms  them  into  a  set  of  first  order  differential  equations  in  single  variables, 
which  are  trivially  simple  to  solve.  The  difficulty  w'ith  this  approach!  is  that, 
because  of  the  z  variation  of  the  elements  of  A(z),  the  diagonalization  equation 
involves  a  "Liapounov  transformation"^^  and  takes  the  form 

(t‘^  at  -  t"^-—).,  =  6.  E..  (39) 

dz  ik  ik  ik 

where  6..  is  the  Kronecker  delta  and  B..  the  i—  element  of  a  matrix, 
ik  ik 

To  find  the  transformation  matrix  T(z)  that  diagonalized  A(z),  it  would  be 
necessary  to  solve  a  set  of  first  order  coupled  ordinary  differential  equations 
that  may  be,  in  the  genera.,  case,  as  complicated  as  the  original  set;  thus  the 
problem  is  not  reduced  to  the  purely  algebraic  one  of  find  the  roots  of  the 
characteristic  equation,  as  it  would  if  A(z)  were  independent  of  z. 
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still  another  analytical  technique  is  the  direct  use  of  the  "matricant  series” 
(11,  12)^  which  evolves  naturally  from  an  iterative  solutijn  of  the  equation 
system  (37)  and  which  is  known  to  converge  to  the  exact  solution.  In  certain 
cases,  e.  g.  where  the  matrix  elements  can  be  represented  by  the  first  tv/o  o’' 
tnree  terms  of  a  power  series  in  z  or  by  a  constant  term  plus  two  or  three 
exponential  or  Fourier  terms,  we  might  obtain  a  series  whose  convergence  is 
fast  enough  to  offer  a  possibility  of  practical  solution.  The  actual  implementa¬ 
tion  of  this  method  to  obtain  useful  solutions  requires  a  high-speed  computer. 

In  the  case  of  self-commutative  matrix  elements  (i.  e.  if  Condition  (38)  above 
holds,)  the  matricant  series  closes,  and,  as  remarked  above,  the  solution  of 
the  equation  system  becomes  tantamount  to  finding  the  eigenvalues  of  the 
matrix  A.  This  point  will  be  discussed  in  further  detail  below. 

The  matricant  method,  although  not  necessarily  practical  in  a  given  case, 
is  nevertheless  worth  discussing  because  of  the  light  it  sheds  on  the  nature  of 
the  solutions.  Consider  the  general  system  of  inhomogeneous  equations  in 
matrix  form 

=  A(z)  u(z)  +  f(z)  (40) 

dz 

where  f(z)  is  a  vector  consisting  of  known  "source”  elements.  Integrating 
(40)  once,  we  have 

u(z)  -  u(z  )  +  dz’  A(z')  u(z’)  +  f(z’)  dz'  .  (41) 

o  z  z 

o  o 

Repeated  substitution  of  u(z)  as  given  by  (41)  into  the  integrand  of  the  first 
integral  on  the  right-hand  side  of  (41)  yields 
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r  z  ^1 ' 

u(z)  =  I  I  +  /^  dz  '  A(z  ')  +  dz  '  A(z  •)  f  dz  '  A(z^') 

L  Zq  1  1  ^  ^  ^ 

Z.  '  z  ' 

+  dz  '  A(z/)  f  ^  dz  '  A(z  ')  /  ^  dz  '  A(z  ')  +  .  .  .  u(z  ) 
o  o  o 

'  -7  ^  ' 

+  4'  +  4^  A(z^')  ^dz^'Hz^') 


+  dz^'  A(z^')  "dz^'  A(z2')  ^dz^'Uz^')  + 


where  I  is  the  identity  matrix. 

Let  us  find  "Green's  function"  of  the  equation  system  (40)  i.e.  the  response 

to  forcing  functions  f(z)  that  are  impulses  at  z  =  z.,  where  z,  >  z  .  We  write 

1  1  o 

f(z)  =  f  6(z  -  z  )  ;  z  <  z  <  z  (43) 

1  1  o  1 

where  f^  is  a  constant,  and  we  choose  at  a  point  where 

u(z  )  -  0  (44) 

o 

which  is  always  allowable  since  the  choice  of  is  perfectly  arbitrary. 


In  this  case,  (4'?)  has  the  form 


z  ' 

u(z)  =  I  +  'iz'  A(z^')  +  dz^'  A(z^')  ^  A(z2')  +  • 

1  11 


.  .  f 

J  X 


From  (42)  and  (45),  we  see  that  solving  the  homogeneous  equation  (37)  for 

a  fixed  set  of  values  of  the  elements  of  u(z  )  is  equivalent  to  solving  the 

o 

inhomogeneous  equation  (40)  with  u(z  )  =  0  and  with  impulse  sources  at  z  =  z  , 
- -  o  1 

where  z^  <  z^  <  z.  The  method  to  be  used  here  is  that  of  assigning  values  to 

u(z^)  and  s  dving  the  homogeneous  equation.  Our  solutions,  then,  will  relate 

to  a  fixed  set  of  values  of  u(z^)  with  the  "forcing  function"  f(z)  set  equal  to 

zero,  and  we  will  dall  z  the  "position  of  the  source". 

o 
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Referring  back  to  (42)  with  f(z)  =  0,  it  is  not  difficult  to  see  that,  if 
the  condition  (38)  holds,  then  the  third  term  of  the  matricant  series  is 


1  .7. 


dz^’A(z^')  dz^'A^lz^')  =  ^  dz^'dz2'A(Zj')A(z2') 


'h(/z 


and  the  general  term  of  the  series  is 

z  '  z'  ^  .  N 

dz^'A(z^');^^  dz2'A(z2')  -"4  '  dZj^A(z^')  =  — dz'A(z')  j  •  (47) 

The  series  then  becomes  th  matrix  analog  of  the  power  series  for  an 
exponential  function  and  is  denoted  by 
A(z')az' 

u(z)  =  (e  °  u(z  !  .  (48) 

V  /  o 

In  the  particular  case  where  the  matrix  elements  are  constant  between 

z  and  z, 
o 

A(z  z  ) 

u(z)  =  e  °  u(z  )  .  (49) 

o 

The  exponential  function  of  a  matrix  B  can  be  evaluated  from  Sylvester's 
Theorem'  '  obtainable  by  Laplace  transform  techniques'  '  and  has  the  ferm 

A(z  -  z  )  n  ^ 

f  /k  '  ^ 

whej  a  the  's  are  the  n  independent  (nor -degenerate)  eigenvalues  of  the  matrix 

iC 

A.  If  degenerate  eigenvalues  exist,  then  a  modification  of  (50)  is  required,  as 

^  1 

discussed  by  Rekoff' 

The  solution  of  (37)  with  condition  (38)  in  effect,  then,  is  given  by  (50), 


with  modifications  if  degenerate  eigenvalues  exist,  and  its  determination 
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requires  that  the  eigenvalues  of  B  be  found.  If  A  is  independent  of  z.  the 
eigenvalues  of  B  are  those  of  A(z-z^);  hence  the  problem  reduces  to  that  of 
finding  the  eigenvalues  of  A.  This  is  the  same  problem  that  must  be  solved 
in  a  dispersion  analysis,  where  the  entire  system  of  equations  is  complex- 
Fourier  or  Laplace  transformed  and  is  therefore  reduced  to  a  homogeneous 
set  of  simultaneous  algebraic  equations.  The  condition  required  for  existence 
of  nontrivial  solutions  of  this  system  is  the  \  mishing  of  the  determinant  of 
coefficients.  Tnis  condition  is  exactly  the  same  as  the  irxdicial  polynomial 
equation,  whose  solution  yields  the  eigenvalues  of  the  A- matrix. 

If  the  elements  of  the  A- matrix  have  certain  simple  z -variations, 
then  solutions  can  be  obtained  through  integral  transforms  or  orthogonal 
function  expa.isions. 

We  have  investigated  a  number  of  such  possibilities  and  found  them 
all  wanting  in  one  way  or  a^iother.  For  example,  a  stratified  medium  theory 
can  be  formulated  in  which  the  z -variation  of  the  elements  of  A(z)  is  linear 
or  exponential.  In  the  former  case,  Laplace  transformation  leads  to  a  first 
order  differential  equation  system  with  coefficients  that  are  linear  in  s.  the 
transform  variable.  Thus  the  problem  is  not  reduced  to  algebra  as  it  is  in 
+he  constant  coeificient  case,  but  rather  to  a  system  of  diffei-ential  equations 
of  the  same  degree  of  difficulty  as  the  original  system.  In  the  case  of 
exponential  elements  of  A(z),  the  transformation  leads  to  difference  equa¬ 
tions  which  can  be  solved  numerically.  Since  inversion  would  be  required 
after  solution,  in  order  to  transform  from  the  s -domain  back  to  the  z -domain, 
this  did  not  seem  like  a  satisfactory  method.  The  same  sort  of  limitations 
seemed  to  exist  in  the  use  of  orthogonal  function  expaiisions.  The  case  of 
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periodic  coefficients,  for  example,  can  be  treated  by  expanding  both  the 
elements  of  A(z)  and  the  solutions  in  Fourier  series  on  z,  but  the  result  is 
a  set  of  coupled  recursion  relations  for  the  coefficients,  instead  of  the 
simultaneous  equations  that  v.'ould  appear  if  the  A-elements  were  constant. 
These  recursion  relations  could  be  solved  but  the  numerical  work  required 
to  obtain  the  Fourier  coefficients  of  the  A(z)  elements  would  be  considerable. 

Another  approach  that  was  attempted  in  the  early  stages  of  investiga¬ 
tion  was  reduction  (by  algebra  and  differentiation)  of  the  original  20-equation 
system  down  to  single  differential  or  integro -differential  equations  of  high 
order  in  single  unknowns.  This  method  iu  feasible  and  perfectly  valid,  but 
it  has  some  very  serious  practical  limitations  if  the  coefficients  have 
arbitrary  z-variation.  Fiist,  the  coefficients  obtained  for  the  resulting 
differential  equations  are  extremely  complicated  algebraically  and  the 
probability  of  eliminating  errors  in  setting  up  compute^’  solutions  would  be 

small.  Secondly,  in  order  to  solve  the  high  order  different’al  equations 
th. 

(e.  g.  12 —  order,  as  in  one  such  development  that  was  carried  out),  it  might 
be  necessary  to  reduce  them  to  sets  of  coupled  first  order  differential  equa¬ 
tions,  v/hich  amounts  to  traveling  in  a  circle. 

Thirdly,  the  coefficients  of  the  equations  finally  obtained  contain  many 
z-derivatives  of  the  ambient  gas  parameters.  Ibis  computation  of  these 
derivatives  would  be  required  before  the  coefficients  of  the  differential  equa¬ 
tions  could  be  specified.  The  derivatives,,  whose  evaluation  would  require 
an  enormous  amount  of  computational  labor,  would  be  highly  inaccurate  since 
they  would  be  obtained  from  very  crude  experimentally  derived  data  curves. 
Thus  a  method  that  contains  as  few  z-deri^'atives  as  possible  in  the  co- 


efficients  is  most  desirable  from  the  viewpoint  of  accurac;;> . 


4.  1.  3.  The  Range -Kutta  Step-by-Step  Method 

No  general  solution  obtainable  by  purely  analytical  methods  was  found 
that  is  suitable  for  the  general  problem,  e.  g.  cases  where  the  z -variation 
of  the  elements  of  A(z)  is  obtained  from  experimental  data  and  where  these 
elements  are  allowed  to  have  any  arbitrary  z -variation  and  do  not  obey  the 
condition  (38).  Such  is  the  case  in  an  accurate  study  of  wave  motions  in  the 
ionosphere;  hence,  we  will  invoke  a  purely  numerical  (stop-by-step)  technique 
which  does  not  lead  to  convergence  problems  even  if  the  elements  have  a 
particularly  complicated  variation  with  z. 

After  studying  the  relative  advantages  of  various  numerical  techniques 
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for  solution  of  the  equation  system'  '  in  the  general  case,  it  was  decided 
to  use  the  Runge-Kutta  step-by-step  procedure.  This  method  has  the 
advantage  of  high  attainable  accuracy  without  a  prohibitively  large  expenditure 
of  computer  time. 

To  apply  the  Runge-Kutta  method  to  (37),  we  first  choose  a  z -interval 
h-  We  then  begin  with  a  set  of  values  of  the  elements  of  the  vector  u(z)  ai 
z  =  z^.  The  i~  element  is  denotea  in  general  by  u^(z).  Its  value  at  z  =  z^  is 

denoted  by  u. 

^  1,  o 

The  long-hand  expression  f  r  the  matrix  equation  (37)  is 


du.(z) 

1 

dz 


N 


L 

j  =  l 


A.  .(z)  u.(z)  ;  i 
J 


=  1, 


(37)' 


z  =  z  +  h.  to 
1  o 


(51) 


■3x 


u  =  u.  +  ^  (k.  +  2k.  +  2k.  ^^^  +  k. 

1,  1  ..  o  6  1,  1  1,  2  1, 3  1, 4 


where 


(1) 


N 


k.  =  h  L  A..(z  )u. 

1.1  IJ  o  J.O 


N 

k.  =  h  Z  A..(z  + 

1.2  13  o  2 

k.  =  h  Z  A..(z  +  5) 

1.3  ij  o  2 


K  1 

u.  +  -iii 
L  hO  2 


(1) 


(1) 


u.  +  ^ 
L  Jif'  2 


(1) 


N 


k.  =  h  Z  A..(z  +h)  u.  +  k. 
1.4  ij  o  Lj.o  3. 


(!)■ 


The  values  of  u^(z)  at  z  =  +  (p  +  l)h,  where  p  is  a  positive  integer,  are 


denoted  by  u.  ,  and  are  given  bv 
1,  p+1 

u.  ^  u.  +  i  (k.  +  2k.  +  2k.  +  k. 

1,  P+1  1,  p  6  1,  1  1,  2  1, 3  1, 4 


(51)' 


where 


k.  =  h  Z  A..(z  +  ph)u. 

1.1  13  o  3,p 


V 

N 

^i  ^  ^  tP  +  i]li)  fu  „  + 

1.  13  o  2  \  3.  P 


(p+1) 


(p+1) 


3.2 


k  =  h 

1,  4 


N 

Z  A..(z  +[p  +  l]h)fu.  +k.  . 

3  =  3  13  o  V  3.P  3.3  J 


Beginning  with  values  of  u.  for  3  =  1,  .  .  .  ,  N,  we  can  proceed  to  find 

3.  ® 

the  successive  values  of  u.  ,,  u.  .  .u,  ,  to  any  desired  number  of 

3.1  32  3, p+1 

steps  by  the  u.'.e  of  eq.  (51)'.  The  accuracy  attainable  by  this  method  is 
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fourth  order  in  h.  The  sharper  the  gradients  in  the  elements  of  A..,  the 
smaller  the  value  of  h  required  for  a  given  level  of  accuracy. 

A  Fortran  program  (soon  to  be  run  on  the  IBM -TOGO)  has  been  prepared 
based  on  the  coupling  of  a  simultaneous  equations  routine  with  a  Runge-Kutta 
routine.  In  effect,  this  program  implements  the  steps  between  eqns.  (35)  and 
(37),  then  solves  (37)  by  the  Runge-Kutta  routine.  We  originally  carried  out 
the  steps  from  (36)  to  (37)  algebraically,  resulting  in  eqn.  (37)  where  the 


^  A, 


elements  of  u(z)  are  the  10  variables  U  ,  U  ,  P  ,,  P  ,,  P, ,,  E  ,  E  ,  H,  ,  J 

nz  pz  nl  pi  il  X  y  lx  z 

and  It  was  then  decided  that  the  direct  computer  solution  of  the  entire 

system  is  more  efficient,  because  some  of  the  11  variables  that  appear  only 
algebraically  (U  ,  U  ,  U  ,  U  >  J  ,  J  .(R  E  ,  B,  )  will  be  desired 

as  outputs,  and  the  computer  can  produce  them  directly  wither difficulty. 

By  means  of  this  program,  it  is  possible  to  study  rnany  types  of  linear  wave 
propagation  in  a  neutral-electron- ion  gas  that  is  non-uniform  in  only  one 
direction.  Any  desired  sets  of  source  functions,  horizontal  propagation  constant 

A 

k  ,  frequency  «  and  z-variation  of  ambient  neutral  and  electron  density, 
collision  frequency  and  static  magnetic  field  magnitude  and  direction  can  be 
studied  by  merely  changing  the  input  numbers.  The  computer  outputs  will 


consist  of  any  desired  components  of  the  velocities  and  U^,  the  current 

A  A  A 

density  J,  the  p'^rturbed  electric  and  magnetic  fields  E  and  B^,  the  densities 

A  A  AAA 

(ft  ,  or  (ft  ,  or  the  pressures  P  ,,  P  ,,  P„ 
nl  pi  ^  nl  pi  il 


4.  2.  Applications 

The  possible  applications  of  the  analysis  discussed  in  this  report  ate 
manifold.  The  computer  program  is  sufficiently  general  to  cover  many  types 
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of  linear  wave  propagation  in  a  neutral  electron- ion  gas.  By  merely  changing 
input  parameter  values,  the  computer  can  produce  answers  to  such  diverse 
questions  as  (a)  the  effect  of  acoustic -gravity  waves  on  properties  of  the 
electron-ion  plasma  in  a  vertically  uihomogeneous  ionosphere,  (b)  the  effects 
of  vertical  non -uniformity  in  electron  density,  collision  frequency  and/or 
ambient  magnetic  field  on  the  propagation  of  MHD  waves  high  in  the  ionosphere 
and  that  of  electromagnetic  waves  in  various  parts  of  the  ionosphere.  In  ".e 
case  of  electromagnetic  waves,  the  point  of  view  to  be  taken  is  that  the  neutral 

A  A  A 

gas  parameters  are  set  to  zero  in  the  system  equations  and 

the  plasma  equations  are  used  to  find  a  conductivity  tensor.  The  latter  is 
vertically  non-uniform,  its  exact  functional  variation  with  z  being  determined 
by  the  z-variation  of  the  static  plasma  parameters.  The  conductivity  tensor 
is  then  substituted  into  the  Maxwell  equations  IV  and  V  in  our  system.  This 
particular  facet  of  vertically  inhomogeneous  ionospheric  wave  +  .eory  has  been 
treated  quite  extensively  in  the  literature.  Examples  of  studies  of  radio  wave 

propagation  with  altitude -varying  conductivity  tensor  are  provided  by  the  work 

(19)  (20)  (21) 

of  Clemmow  and  Heading  ,  Budden  and  Clemmow  ,  Barron  and  others. 

The  acoustic -gravity  wave  plasma  interaction  problem  is  based  on  an 
entirely  different  point-of-view.  Over  most  of  the  ionosphere,  the  fractional 
ionization  is  so  small  that  the  terms  in  the  neutral  gas  equations  relating  to  the 
effects  of  the  plasma  are  negligible,  except  at  frequencies  that  are  extremely 
low  relative  to  neutral-plasma  collision  frequencies.  The  neutral  gas  wave 

This  reflects  the  fact  that,  because  of  low  fractional  ionization  in  the 
ionosphere,  an  electromagnetic  wave  would  not  significantly  perturb  the 
neutral  gas. 


prop -’’.gallon  is  thus  only  very  weakly  influenced  by  the  presence  of  plasma. 

The  neutral  waves, whose  equations  are  now  uncoupled  from  the  plasma  equa- 

(3) 

tions,  can  by  analyzed  by  the  methods  of  Hines'  19v30  paper  and  the  result- 

^  A. 

ing  neutral  gas  parameters  U  .  P  ,  and  (ft  ,  can  then  be  used  as  source  terms 

j^nl  nl  nl 

in  the  plasma  equations.  In  this  way,  the  system,  which  in  its  general  form  is  a 
system  of  10  homogeneous  differential  equations  and  11  homogeneous  algebraic 
equations  in  a  total  of  21  unknowns,  now  becomes  a  system  of  8  inhomogeneous 
differential  equations  and  8  algebraic  equations  in  a  total  of  16  unknowns  and 
with  known  source  terms  determined  by  solving  the  first  5  equations  of  the 
original  system. 

This  technique,  which  takes  advantage  of  the  extremely  low  value  of  a 
that  prevails  in  the  ionospnere,  reduces  the  magnitude  of  our  computer  problem 
in  treating  this  particular  ionospheric  effect.  A  detailed  discussion  of  the 
acoustic  gravity  wave  plasma  interaction  problem  will  be  presented  in  a  later 
report. 

Our  computer  program  is  by  no  means  restricted  to  the  ionosphere.  It 
can  be  applied  to  any  linear  neutral-electron-ion  gas  with  static  parameters 
that  are  non-uniform  in  a  single  direction.  The  normalization  of  the  parameters 
in  our  equations  wou'.d  enhance  the  convenience  of  studying  such  gases,  i.  e.  , 
the  normalized  parame':ers  could  be  varied  in  such  a  way  as  to  produce 
universal  curves.  Wave  problems  that  involve  complicated  one -dimensional 
spatial  variations  of  parameters,  simple  degenerate  cases  of  which  have  been 
solved  by  purely  analytical  methods,  can  be  treated  by  our  computer  program 
without  the  necessity  for  either  drastically  simplifying  the  parameter  varia¬ 
tions  or  resorting  to  perturbation  techniques.  The  latter,  of  course,  are 

Provided,  of  course,  that  the  assumptions  and  approximations  used  here 
are  applicable. 
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severely  limited  in  accuracy  when  the  departure  from  the  idealized  case  is 
large . 
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TABLE  1 


Ordering  of  Equations  and  Vectors  in  Equations  (34)  and  (35) 
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u. 
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TABLE  2 

Matrices  in  Equations  34  and  35 
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TABLE  2  (cont'd) 


Matrices  in  Equations  34  and  35 
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